
divlimit <- 0.3
name <- "sim" 

pdf("divergence.632.pdf",width=6,height=6,pointsize=12)
par(mar=c(3,3,3,1))
for(i in 632:632){

      	params <- read.csv(paste(name,i,".in",sep=""),header=FALSE)
	par(pty="s")
	divergence <- read.csv(paste(name,i,"/",name,i,".divergence.log",sep=""),header=FALSE)
	nmu <- format(1/as.numeric(levels(params[22,1]))[params[22,1]],scientific=TRUE)
	smu <- format(1/as.numeric(levels(params[23,1]))[params[23,1]],scientific=TRUE)
	
	plot(divergence[,1]/364,divergence[,2],xlab=NA,ylab=NA,t="l",ylim=c(0,divlimit),xlim=c(0,7300/364),col="black")
	fixations <- read.csv(paste(name,i,"/",name,i,".fixations.log",sep=""),header=FALSE)
	fixations <- cbind(fixations[,1:4], divergence[,2])
	# neutral fixations by themselves, or together with repro and survival	
	fixev <- fixations[fixations[,2]>fixations[,3]+fixations[,4],]
	points(fixev[,1]/364,fixev[,5],col="black",pch="|")
	# survival fixations
	fixev <- fixations[fixations[,3]>0,]	
	points(fixev[,1]/364,fixev[,5],col="black",pch="X")
	# reprod fixations
	fixev <- fixations[fixations[,4]>0,]
	points(fixev[,1]/364,fixev[,5],col="black",pch="O")
	# Write number of clones that went to fixations at that point in time
	if(nrow(fixev)>0){text(fixev[,1]/364,fixev[,5]+0.02,fixev[,4])}
	# generations
	f <- divergence[nrow(divergence)/4*c(1:5),3]
	axis(side=3,at=20*c(1:5)/4,labels=f)
	rval <- t(params)[31]
	sval <- t(params)[15]
	mtext(bquote(theta == (.(i) ~ "," ~ .(rval) ~ "," ~ .(sval) ~ "," ~ .(nmu) ~ "," ~ .(smu))), side=1, line = 2)
	mtext("Divergence", side=2, line = 2)
	mtext("Generations", side=3, line = 2)
}

dev.off()





